simulation results/prelim_gamma_results.R

In BriceonWiley/IntegratedLikelihood.R:

library(parallel)

gam_50_250 <- readRDS("~/Documents/Coding/RStudio/Research/IntegratedLikelihood.R/simulation results/gam_50_250.rds")
gam_300_500 <- readRDS("~/Documents/Coding/RStudio/Research/IntegratedLikelihood.R/simulation results/gam_300_500.rds")

# for plotting
DATA %>%
  map_dfr(~ coverage_gam(.x, TRUE)) %>%
  select(-N) %>%
  gather("Interval", "Coverage", -gamma, -p1, -p2) %>%
  ggplot(aes(gamma, Coverage)) +
  geom_path(aes(color = Interval)) +
  geom_hline(yintercept = 0.95, linetype = "dashed") +
  theme_minimal()

# simulation 50 to 100
mix <- expand_grid(p = seq(0.05, 0.95, by = 0.02), n = seq(50, 100, by = 50))
# sim_big_50_100 <- mcMap(
#   function(N, p) simulate_gamma(
#     N1 = N, p1 = p, phi = 0.1, theta1 = 0.15,
#     p2 = 0.3, phi2 = 0.2, theta2 = 0.25, len = 10000
#   ),
#   N = mix$n, p = mix$p, mc.cores = 4
# )
# begin 3/21 @ 9:00
sim_big_250_c <- mclapply(
  seq(0.05, 0.95, by = 0.05),
  function(p) simulate_gamma(
    N1 = 250, p1 = p, phi1 = 0.1, theta1 = 0.15,
    p2 = 0.3, phi2 = 0.2, theta2 = 0.25, len = 2500
  ),
  mc.cores = 4
)
# for saving on desktop
saveRDS(sim_big_250_c, "/Users/briceon_wiley/Documents/RStudio/IntegratedLikelihood.R/simulation results/results_gam/gam_250_big_part_3.rds")

sim_plots <-
  map2(
    1:5, seq(50, 250, by = 50),
    ~ {
      sim_full %>%
        map_dfr(~ coverage_gam(.x, TRUE)) %>%
        group_split(N) %>%
        pluck(.x) %>%
        gather("Interval", "Coverage", -gamma, -p1, -p2, -N, -NaNs) %>%
        ggplot(aes(gamma, Coverage)) +
        geom_path(aes(color = Interval)) +
        geom_hline(yintercept = 0.95, linetype = "dashed") +
        theme_minimal() +
        ylim(c(0.5, 1)) +
        labs(
          title = glue("N = {.y}"),
          subtitle = "theta1 = 0.15, p2 = 0.30, phi2 = 0.20, theta2 = 0.25"
        ) +
        geom_text(
          aes(
            x = 0, y = 0.6,
            label = glue(
              "ILR Coverage: min = {round(min(Coverage[Interval == 'ILR']), 4)}, max = {round(max(Coverage[Interval == 'ILR']), 4)}"
            )
          )
        ) +
        geom_text(
          aes(
            x = 0, y = 0.55,
            label = glue(
              "W Coverage: min = {round(min(Coverage[Interval == 'W']), 4)}, max = {round(max(Coverage[Interval == 'W']), 4)}"
            )
          )
        )
    }
  )

map2(
  sim_plots, seq(50, 250, by = 50),
  ~ ggsave(glue("gam_plot_{.y}.pdf"), .x)
)